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Abstract 

In this contribution, we develop a variational integrator for the simulation of (stochastic and multiscale) electric 
circuits. When considering the dynamics of an electrical circuit, one is faced with three special situations: 1. The 
system involves external (control) forcing through external (controlled) voltage sources and resistors. 2. The system 
is constrained via the Kirchhoff current (KCL) and voltage laws (KVL). 3. The Lagrangian is degenerate. Based on 
a geometric setting, an appropriate variational formulation is presented to model the circuit from which the equations 
of motion are derived. A time-discrete variational formulation provides an iteration scheme for the simulation of 
the electric circuit. Dependent on the discretization, the intrinsic degeneracy of the system can be canceled for the 
discrete variational scheme. In this way, a variational integrator is constructed that gains several advantages compared 
to standard integration tools for circuits; in particular, a comparison to BDF methods (which are usually the method 
of choice for the simulation of electric circuits) shows that even for simple LCR circuits, a better energy behavior and 
frequency spectrum preservation can be observed using the developed variational integrator. 

Keywords: structure-preserving integration, variational integrators, degenerate systems, electric circuits, noisy 
systems, multiscale integration 



1. Introduction 

Variational integrators have mainly been developed and used for a wide variety of mechanical systems. However, 
real-life systems are generally not of purely mechanical character. In fact, more and more systems become multidis- 
ciplinary in the sense, that not only mechanical parts, but also electric and software subsystems are involved, resulting 
in mechatronic systems. Since the integration of these systems with a unified simulation tool is desirable, the aim of 
this work is to extend the applicability of variational integrators to mechatronic systems. In particular, as the first step 
towards a unified simulation, we develop a variational integrator for the simulation of electric circuits. 

Overview. Variational integrators [30] are based on a discrete variational formulation of the underlying system, e.g. 
based on a discrete version of Hamilton's principle for conservative mechanical systems. The resulting integrators 
given by the discrete Euler-Lagrange equations are symplectic and momentum-preserving and have an excellent long- 
time energy behavior. Choosing different variational formulations (e.g. Hamilton, Lagrange-d'Alembert, Hamilton- 
Pontryagin, etc.), variational integrators have been developed for classical conservative mechanical systems (for an 
overview see [24, 25]), forced [18] and controlled [34] systems, constrained systems (holonomic [26, 27] and non- 
holonomic systems [19]), nonsmooth systems [11], stochastic systems [6], and multiscale systems [39]. Most of 
these systems share the assumption, that they are non-degenerate, i.e. the Legendre transformation of the correspond- 
ing Lagrangian is a diffeomorphism. Applying Hamilton's principle to a regular Lagrangian system, the resulting 
Euler-Lagrange equations are ordinary differential equations of second order and equivalent to Hamilton's equations. 
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The Lagrangian formulation for LC circuits is based on the electric and magnetic energies in the circuit and the 
interconnection constraints expressed in the Kirchhoff laws. There exist a large variety of different approaches for a 
Lagrangian or Hamiltonian formulation of electric circuits (see e.g. [9, 21, 10, 2, 38] and references therein). All of 
theses authors treat the question of which choice of the Lagrangian coordinates and derivatives is the most appropriate 
one. Several settings have been proposed and analyzed, e.g. a variational formulation based on capacitor charges 
and currents, on inductor fluxes and voltages, and a combination of both settings, as well as formulations based on 
linear combinations of the charges and flux linkages. Typically, one wants to find a set of generalized coordinates, 
such that the resulting Lagrangian is non-degenerate. However, within such a formulation, the variables are not easily 
interpretable in terms of original terms of a circuit. 

A recently-considered alternative formulation is based on a redundant set of coordinates resulting in a Lagrangian 
system for which the Lagrangian is degenerate. For a degenerate Lagrangian system, i.e. the Legendre transform is 
not invertible, the Euler-Lagrange equations involve additional hidden algebraic constraints. Then, the equations do 
not have a unique solution, and additional constraints are required for unique solvability of the system. For the circuit 
case, these are provided by the Kirchhoff Current Law (KCL). From a geometric point of view, the KCL provides a 
constraint distribution that induces a Dirac structure for the degenerate system. The associated system is then denoted 
by an implicit Lagrangian system. In [42] and [41], it was shown that nonholonomic mechanical systems and LC 
circuits as degenerate Lagrangian systems can be formulated in the context of induced Dirac structures and associated 
implicit Lagrangian systems. The variational structure of an implicit Lagrange system is given in the context of the 
Hamiltonian-Pontryagin-d'Alembert principle, as shown in [43]. The resulting Euler-Lagrange equations are called 
the implicit Euler-Lagrange equations [42, 43, 32], which are semi-explicit differential-algebraic equations that con- 
sist of a system of first order differential equations and an additional algebraic equation that constrains the image of the 
Legendre transformation (called the set of primary constraints). Thus, the modeling of electric circuits involves both 
primary constraints as well as constraints coming from Kirchhoff' s laws. In [17], an extension towards the intercon- 
nection of implicit Lagrange systems for electric circuits is demonstrated. For completeness, we have to mention that 
the corresponding notion of implicit Hamiltonian systems and implicit Hamiltonian equations was developed earlier 
by [3, 31, 4]. An intrinsic Hamiltonian formulation of dynamics of LC circuits as well as interconnections of Dirac 
structures have been developed, e.g. in [31] and [8], respectively. 

There are only a few works dealing with the variational simulation of degenerate systems, e.g. in [36], variational 
integrators with application to point vertices as a special case of degenerate Lagrangian system are developed. Al- 
though there exists a variety of different variational formulations for electric circuits, variational integrators for their 
simulation have not been concretely investigated and applied thus far. In [22], a framework for the description of the 
discrete analogues of implicit Lagrangian and Hamiltonian systems is proposed. This framework is the foundation 
for the development of an integration scheme. However, no concrete simulation scenarios have yet been performed. 
Furthermore, the discrete formulation of the variational principle is slightly different from the approach presented in 
this work, thus resulting in a different scheme. 

Contribution. In this work, we present a unified variational framework for the modeling and simulation of electric 
circuits. The focus of our analysis is on the case of ideal linear circuit elements, consisting of inductors, capacitors, 
resistors and voltage sources. However, this is not a restriction of this approach, and the variational integrators can also 
be developed for nonlinear circuits, which is left for future work. A geometric formulation of the different possible 
state spaces for a circuit model is introduced. This geometric view point forms the basis for a variational formula- 
tion. Rather than dealing with Dirac structures, we work directly with the corresponding variational principle, where 
we follow the approach introduced in [43]. When considering the dynamics of an electric circuit, one is faced with 
three specific situations that lead to a special treatment within the variational formulation and thus the construction 
of appropriate variational integrators: 1. The system involves external (control) forcing through external (controlled) 
voltage sources. 2. The system is constrained via the Kirchhoff current (KCL) and voltage laws (KVL). 3. The La- 
grangian is degenerate leading to primary constraints. For the treatment of forced systems, the Lagrange-d' Alembert 
principle is the principle of choice. Involving constraints, one has to consider constrained variations resulting in a 
constrained principle. The degeneracy requires the use of the Pontryagin version; thus, the principle of choice is 
the constrained Lagrange-d' Alembert-Pontry agin principle [43]. Two variational formulations are considered: First, 
a constrained variational formulation is introduced for which the KCL constraints are explicitly given as algebraic 
constraints, whereas the KVL are given by the resulting Euler-Lagrange equations. Second, an equivalent reduced 
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constrained variational principle is developed, for which the KCL constraints are eliminated due to a representation of 
the Lagrangian on a reduced space. In this setting, the charges and flux linkages are the differential variables, whereas 
the currents play the role of algebraic variables. The number of inductors in the circuit and the circuit topology de- 
termine the degree of degeneracy of the system. For the reduced version, we show for which cases the degeneracy 
of the system is canceled via the KCL constraints. Based on the variational formulation, a variational integrator for 
electric circuits can be constructed. For the case of a degenerate system, the applicability of the variational integrator 
is dependent on the choice of discretization. Based on the type and order of the discretization, the degeneracy of the 
continuous system is canceled for the resulting discrete scheme. Three different integrators and their applicability to 
different electric circuits are investigated. The generality of a unified geometric (and discrete) variational formulation 
is advantageous for the analysis - for very complex circuits in particular. Using the geometric approach, the main 
structure-preserving properties of the (discrete) Lagrangian system can be derived. In particular, good energy behav- 
ior and preservation of the spectrum of high frequencies of the solutions can be observed. Furthermore, preserved 
momentum maps due to symmetries of the Lagrangian system can be derived. Going one step further, we extend the 
approach to a stochastic and multiscale setting. Due to the variational framework, the resulting stochastic integrator 
will well capture the statistics of the solution (see for instance [7]), and the resulting multiscale integrator will still be 
variational [39]. 

Outline. In Section 2, we first review the basic notation for electric circuits followed by a graph representation to 
describe the circuit topology. In addition, we introduce a geometric formulation that gives an interpretation of the dif- 
ferent state spaces of a circuit model. Based on the geometric view point, the two (reduced and unreduced) variational 
formulations are derived in Section 3. The equivalence of both formulations as well as conditions for obtaining a non- 
degenerate reduced system are proven. In Section 4, the construction of different variational integrators for electric 
circuits is described and conditions for their applicability are derived. The main structure-preserving properties of the 
Lagrangian system and the variational integrator are summarized in Section 5. In Section 6, the approach is extended 
for the treatment of noisy circuits. In Section 7, the efficiency of the developed variational integrators is demonstrated 
by means of numerical examples. A comparison with standard circuit modeling and circuit integrators is given. In 
particular, the applicability of the multiscale method FLAVOR [39] is demonstrated for a circuit with different time 
scales. 

2. Electric circuits 

2.1. Basic notations 

Considering an electric circuit, we introduce the following notations (following [33]): A node is a point in the 
circuit where two or more elements meet. A path is a trace of adjacent elements, with no elements included more than 
once. A branch is a path that connects two nodes. A loop is a path that begins and ends at the same node. A mesh 
(also called fundamental loop) is a loop that does not enclose any other loops. A planar circuit is a circuit that can be 
drawn on a plane without crossing branches. 

Let q(t), v(f), u(t) e W be the time-dependent charges, the currents and voltages of the circuit elements with 
t € [0,T], where qj(t),Vj(t),Uj(t) € W J , J e \L,C,R,V} are the corresponding quantities through the n L inductors, 
the n c capacitors, the n R resistors, and the n v voltage sources. In addition, we give each of those devices an assumed 
current flow direction. In Table 1, the characteristic equations for basic elements are listed. For the analysis in this 
work, we focus on ideal linear circuit elements, resulting in the following constitutive laws: 

u L (t) = Lv L (t), v c (t) = Cu c (f), u R (t) = Rv R (t), 

with inductance L, capacitance C and resistance R = G l with conductance G and where in general it holds q(t) = v(t). 
The flux linkage for each element is denoted by p(t) e K" and for an inductor, it is defined as the time integral of the 
voltage across the inductor. Note that in the case of an inductor (resp. a capacitor), the associated charge q L (resp. flux 
linkage p c ) is an artificial variable. Similarly, for the resistors and the voltage sources, the associated charges q R , q v 
and flux linkages p R ,pv are artificial variables. 
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Table 1: Characteristic equations for basic circuit elements. 



Ideal inductors and capacitors are purely reactive, i.e. they dissipate no energy. Thus, the magnetic energy stored 
in one inductor with inductance L is 



1 2 
F — —Tv 



where the amount of energy storage in one capacitor with capacitance C is 



u c dq= - dq = - -q 2 c 
o=0 Jo=0 <- L <- 



2.2. Graph representation 

Consider now a circuit as a connected, directed graph with n edges and m + 1 nodes. On the ;'th edge, there are: a 
capacitor with capacitance C ; , an inductor with inductance L,, a voltage source e, and a resistor with resistance R h one 
or several of which can be zeros (cf. Figure 1). Thus, branches in the circuit correspond to edges in the graph. In the 
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Figure 1: A typical branch of a circuit. On this edge, there are: an inductor L t , a capacitor C,, a resistor Rj, and a voltage source e,-, one or several 
of which can be zeros. 

special case that each edge in the graph represents only one circuit element, the number of edges in the graph equals 
the number of circuit elements, and the number of nodes of the circuit and the graph are the same. For simplicity, we 
use the notions from circuit theory, i.e. talking about branches and meshes in the graph. 
For the analysis with circuits, one is faced with the following two basic laws: 

1 . The Kirchhoff Current Law (KCL) states that the sum of currents leading to and leaving from any node is equal 
to zero. 

2. The Kirchhoff Voltage Law (KVL) states that the sum of voltages along each mesh (or fundamental loop) of the 
network is equal to zero. 

Let K G R"' m be the Kirchhoff Constraint matrix of a given circuit represented via a graph defined by 

!-l branch i connected inward to node j 
+ 1 branch ;' connected outward to node j (1) 
otherwise. 

In the special case where the two ends of an edge are connected to the same node, we set K t j = 0. Since the ground 
node is excluded, the Kirchhoff Constraint matrix has only m rather than m + 1 columns. Allowing only one circuit 
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element for one branch, either inductor, capacitor, resistor, or voltage source, K can be expressed as 



K = 



K L 
K c 
K R 
K v 



where Kj e R" J ' m , J e \L, C,R, V} is the Constraint Matrix for the set of n L inductors, n c capacitors, n R resistors, and 
riy voltage sources, respectively with n L + n c + n R + n v = n. The Kirchhoff Constraint Matrix provides the Kirchhoff 
current constraints as K T v = 0. For connected, planar graphs, the number of meshes I is determined via I = n - m, 
where n is the number of branches and m + 1 the number of nodes. This is a direct consequence from Euler's formula 
[13]. We can thus define the Fundamental Loop matrix K 2 e M. n '"~ m by 



K 2 jj - 



where again K 2 can be expressed as 



-1 branch i is a backward branch in mesh j 
+ 1 branch i is a forward branch in mesh j 
branch i does not belong to mesh j, 



(2) 



K 2 = 



K 2 ,l 
K 2 ,c 

K 2 ,R 

K 2 v 



with K 2 j e M. nj,n ~' n , J e [L, C, R, V] is the Loop Matrix for the set of Ul inductors, nc capacitors, n R resistors, and ny 
voltage sources, respectively. The Fundamental Loop Matrix provides the Kirchhoff voltage constraints as K% u = 0. 
An alternative expression of the Kirchhoff voltage constraints is given by Kh = u, where u are the node voltages of 
the circuit. By u e ker(Kl) and u e im{K) it follows directly ker(Kl) = im(K) and thus im(K 2 ) ± im(K). 



2.3. Geometric setting 

Using a geometric approach for analyzing circuits, we define the configuration manifold to be the charge space 
Q c W of circuit branches with points on the manifold denoted by q e Q. For a particular charge configuration q, the 
tangent bundle TQ is the current space with currents v e T q Q c M." passing through the branches. The corresponding 
cotangent bundle T*Q is the flux linkage space with the flux linkages p e T*Q c E". Note that due to the analogon 
of the quantities, configuration, velocity and momentum in mechanical systems we stick with the notation (q, v, p) for 
charge, current and flux linkage. The branch voltages u are the analogon of forces for the mechanical system and are 
thus assumed to be covectors in the cotangent space T*Q. 

Let Ag c TQ be a constraint distribution, which is locally given by 

A Q (q) = {v e T q Q I {w a , v) = 0,a=l,...,m}cT q Q (3) 

with the natural pairing (■, ■) : T*Qx T q Q — > R of cotangent and tangent vectors, w" are m independent one-forms 
that form the basis for the annihilator A° Q (q) <zT*Q which is locally given by 

A° (?) = {w e r q Q I (w, v) = Vv e A Q (q)} c T q Q. (4) 

Using the matrix K T as local coordinate representation for the one-forms w", the distribution (3) forms the con- 
straint KCL space given by the submanifold 

A Q (q) = {veT q Q\K T v = 0}cT q Q, 

that is spanned by ker(K T ). Its annihilator A^(q) can thus locally expressed by the image im(K) of K. Chooing this 
coordinate representation and with ker(K^) = im(K), the annihilator (4) decribes the constraint KVL space by 

A° Q (q) = \u er q Q\K T 2 u = 0} <zT;Q 
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Note, that the choices of K and K 2 are in general not unique. The only design criterium for K 2 is the condition 
im(K) ± im(K 2 ). Alternative to (2), a matrix K 2 can be constructed using a QR-decomposition of K. Thus, this 
approach is not restricted to cases, where the mesh topology is obvious as for planar graphs. However, in the following 
we work with the Fundamental Loop matrix as candidate for the matrix K 2 due to the physical interpretation. 

From a geometric point of view we can distinguish between three different spaces: Let 03 denote the space of 
branches, 9Jt the space of meshes and 91 the space of nodes, where we exclude the one node defined as ground. 
(q,v,p,u) denote the branch charges, currents, flux linkages, and voltages, and (q,v,p,u) and (q,v,p,u) the corre- 
sponding quantities in mesh and node space, respectively. From KCL and KVL we know that the node currents (and 
charges) as well as the mesh voltages are zero. For 371 and 9T, we define the corresponding configuration, tangent and 
cotangent spaces M c R"-' n , T q M c R"-' n , T*M c R"-' n and N c R m , T q N c R m , T*N c R m . Then, branch, loop 
and node space are defined to be the Pontryagin bundle consisting of the direct sum of tangent and cotangent space, 
i.e. 03 = A Q © A ( ^, 9Jt = TM © T*M and 9T = TN © T*N. 

The following diagram gives the relation between the defined spaces in terms of the Kirchhoff Constraint matrix 
K and the Fundamental Loop matrix K 2 

T*N — A° — TM 

nodes branches meshes (5) 

TN < An < TM 

K T K 2 



with the linear maps K T : A Q -» TN and K 2 : TM -» A e , and their adjoints A" : T*A^ -> A° and K\ : A° Q -» T*M. As 
already stated above, the branch currents consistent with KCL are determined by ker{K T ), where the branch voltages 
consistent with KVL are given by ker(K^). On the other hand, from diagram (5) we can directly follow, that the set 
of branch currents fulfilling the KCL can alternatively expressed as v = K 2 v, whereas the set of branch voltages that 
fulfill the KVL are in the image of K as u = Ku. These are the standard relations between branch currents v(f) and 
mesh currents v(f) and branch voltages u(t) and node voltages u{t), respectively, given by KCL and KVL. 

Note, that diagram (5) represents the general relations between the tangent bundles TQ, TN, and TM, and the cor- 
responding cotangent bundles. The matrices K and K 2 are local coordinate choices for the embeddings of the different 
spaces. These are not unique, however using different coordinates representations, the submanifolds TN,T*N,TM, 
and T*M loose their physical meaning. 

Following the lines of [43] the tangent space at q can be splitted such that T q Q = "H q © <V q , where < H q = A q (q) is 
the horizontal space and *V 9 the vertical space at q. The matrix K T is a local matrix representation of the Ehresmann 
connection A q :T q Q^> <V q . 

Remark 1. A branch can consist of more than one circuit element in a row. In this case, the branch voltage is assumed 
to be the sum of the voltages of all elements in this branch. 

3. Variational formulation for electric circuits 

In the following, we derive the equations of motion for the circuit system, making use of variational principles 
known in mechanics. Due to the constrained system (via KCL and KVL), we present two different variational formu- 
lations that distinguish the way the constraints are involved. 

3.1. Constrained variational formulation 

We can define a Lagrangian X : TQ — > R of the circuit system consisting of the difference between magnetic and 
electric energy as 

£(q,v)= l -v T Lv- l -q T Cq (6) 
6 



with L = diag(Li, . . . , L n ) and C = diag . . . , In the case where no inductor (resp. no capacitor) is on branch i, 
the corresponding entry L, (resp. ^) in the matrix L (resp. C) is zero. In the presence of mutual inductors rather than 
self inductors, the matrix L is not diagonal anymore, but always positive semi -definite. If not explicitly mentioned the 
following theory and construction is also valid for mutual inductors. The Legendre transform FX : TQ — > T*Q is 
defined by 

¥£(q,v) = (q,d£/dv) = (q,Lv). (7) 

Note that the Lagrangian can be degenerate if the Legendre transform is not invertible, i.e. L is singular. The constraint 
flux linkage subspace 1 is defined by the Legendre transform as 

P = F£(A e )c T*Q, 

where A e c TQ is the distribution. The Lagrangian force of the system consists of a damping force that results from 
the resistors and an external force being the voltage sources 

f L (q, v, f) = -dmg(R)v + diag(£) M (8) 

with /? = (/?!,..., R n ) T and £ = (e\, . . . , e n ) T . If no resistor is on branch ;, the corresponding entry R t in the vector R is 
zero. For the entries of the vector S, it holds e, = if no voltage source is on branch i and e, = 1 otherwise. Here we 
assume that the time evolution of the voltage sources is given as time dependent function u s (t). Thus, in the following, 
we replace diag(£)« by u s (f) for a given function u s : [0, T] — > K". 

To derive the equations of motion for the circuit system, we make use of the Lagrange-d'Alembert-Pontryagin 
principle, i.e. we are searching for curves q(t), v(t) and p(t) fulfilling 

6 f £{q{i), v(0) + (Pit), q(t) - v(0> dt + f f L (q(t\ K0, ' *9(0 * = (9) 
Jo Jo 

with fixed initial and final variations 6q(0) = 6q(T) = and constrained variations 5q € Agiq). 
Taking variations gives us 



a 



+ h, Sq^j - (p, 5q) + (dp, q-v) + (j^ | - p. Sv 



dt = (10) 



for arbitrary variations 5v and dp, K T v = and constrained variations 5q G Ag(<7). This leads to the constrained 
Euler-Lagrange equations 

^-p + f L eA° Q (q) (11a) 

q = v (lib) 

^-P = (11c) 

K T v = 0. (lid) 

For the Lagrangian (6) and the forces (8), the constrained Euler-Lagrange equations are 

p = -Cq - diag(/?)v + u s + KA (12a) 

q = v (12b) 

p - Lv (12c) 

K T v = 0, (12d) 



also denoted by the set of primary constraints 
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where A represent the node voltages u e T*N. Thus, the first line corresponds to the KVL equations of the form 
Ku = u, and the last line are the KCL equations. System (12) is a differential-algebraic system with differential 
variables q and p and algebraic variables v and A. The involvement of the function u s (t) makes the system a non- 
autonomous system. Equation (12c) (also denoted by primary constraints) reflects the degeneracy of the Lagrangian 
system: since FX is not invertible (i.e. L is singular), we can not eliminate the algebraic variable v to obtain a purely 
Hamiltonian formulation. However, in the next step, we eliminate the algebraic variable A by the use of a reduced 
constrained variational principle. 

3.2. Reduced constrained variational formulation 

With the following reduced principle, we derive a slightly different form of the resulting differential-algebraic 
system. This reduced formulation is advantageous from different perspectives: First, the reduced formulation is less 
redundant, such that the Lagrange multipliers are eliminated and the state space dimension is reduced. Second, for 
specific circuits, the degeneracy of the Lagrangian is canceled. Third, the reduced state space still has a physical and 
geometric interpretation: The reduced Lagrangian is defined on the mesh space TM c R 2 ("- m ) rather than on the 
branch space TQ c M. 2 ". 

For the reduction, instead of treating the KCL as extra constraint in the form K T v = 0, we directly involve the 
KCL form K 2 \> = v with v e T q M c W. n ~ m for the definition of the new Lagrangian system. Since K is constant, the 
constraints are integrable, i.e. the configurations q are constrained to be in the submanifold 

C = {qeQ\K T q = 0) 

for consistent initial values qo e C. This simply means that topological relationships that apply for currents will also 
hold for charges to within a constant vector. Then, it holds T q C = A Q (q) and the branch charges q can be expressed 
by the mesh charges q e M c W~ m as q = K 2 q. We define the constrained Lagrangian £ M : TM -» R via pullback 
as t M := K*£ : TM — > R with 

£ M (q, v) = £(K 2 q, K 2 v) = l -v T K T 2 LK 2 v - l -q T K T 2 CK 2 q (13) 

with the Legendre transform ¥£ M : TM -> T*M 

¥£ M (q,v) = (g,d£ M /8v) = (q,KlLK 2 v). 

Depending on the inductor matrix L and the graph topology, the matrix K 2 LK 2 can still be singular, i.e. the Lagrangian 
system can still be degenerate. The cotangent bundle T*M is given by 

TM = {(q, p) e M"- m - n - m l (q, p) = F£ M (q, v) with (q, v) e TM} 
= {(q, p) e R"-™."-" 1 1 (q, p) = (q, K\p) with p e P). 

Thus, the constrained force in T*M is defined as 

fif (q, v, t) = K T 2 f L {K 2 q, K 2 v, t) = -K T 2 A^g{,R)K 2 v + K T 2 u s {t). (14) 

With p e T* q M c R n ~ m given as p = K 2 p we obtain the following reduced Lagrange-d' Alembert-Pontryagin principle 

6 f £ M {q{t), Ht)) + {p(t), ~qit) ~ Kt)) dt + f ffiqit), v(f), f) ■ 6~q{f) dt = (15) 
with fixed initial and final variations 6q(0) = 6q(T) = 0. Taking variations gives us 

X [(^ + ^'^)~ M) + ( sp '^~ 9 ) + ((^w)~ p ' S9 ) dt = (16) 
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for arbitrary variations dv and 5p and 6q. This results in the reduced Euler-Lagrange equations 



d£ M 



P + /l=0 (17a) 



dq 

q = v (17b) 

d£ M 

~w~ p = °- (17c) 

For the Lagrangian (13) and the forces (14), the constrained Euler-Lagrange equations are 

p = Kl (-CK 2 q - diag(R)K 2 v + u s ) (18a) 
q = v (18b) 
p = KlLK 2 v. (18c) 

Here, the first equation is now the KVL in the form K^u = 0, in which the KCL in the form K 2 v — v is also 
involved. System (18) is a differential-algebraic system with differential variables q and p and algebraic variables v. 
The algebraic equation (18c) is the Legendre transformation of the system. If this is invertible (i.e. the matrix K^LK 2 
is regular), the algebraic variable v can be eliminated. In this case, the Euler-Lagrange equations (18) represent a 
non-degenerate Lagrangian system. 

In the following proposition we show for which cases the reduced Lagrangian system is non-degenerate for LC 
circuits, i.e. for which cases the KVL cancels the degeneracy. The statements for RCL and RCLV circuits can be 
derived in an analogous way (see Remark 2b)). 

Proposition 1. For LC circuits (including only self inductors), the system is non-degenerate if the number of capaci- 
tors equals the number of independent constraints involving the currents through the capacities' branches. 

Proof. We have to show that ker(K^ LK 2 ) = {0}. Let «c be the number of capacitors, m the number of Kirchhoff 
Constraints such that K T C e M.' n,nc . Let Iq < m be the number of independent constraints involving the currents through 

With *. /_ ,T70 howa rnwMV^ 

ker{L) = {veT q Q\v L = 0) 



the capacitives branches. With « c = / c < m we have rank(Kl) = n c , thus ker(K^,) = {0}. On the other hand, it holds 



and 

<R{K 2 ) = ker(K T ) = [v e T q Q \ K T v = 0} = [v e T q Q \ (K T L K T C ) | ^ j = 0} 
With ker(K T c ) = {0} this results in 

K(K 2 ) n ker(L) = {v e T q Q \ K T c v c = 0} = {0}. (19) 

and thus ker(LK 2 ) = {0}. Since L is a diagonal matrix, we can split K^LK 2 into yfl? y[LK 2 , where y[Z corresponds 
to the diagonal matrix with diagonal elements VL*. Lj > 0, i = 1, . . . , n. Since K 2 has full column rank, we know with 
(19) that y[LK 2 also has full column rank. It follows for y e R"-' n and y e ker(K^LK 2 ) 

K^LK 2 y = => y T KlLK 2 y = <^ y T K\ VZ r ^LK 2 y = <=> || ^LK 2 y\\ 2 = 

and thus y = since ker{ yfLK 2 ) = {0}. We therefore have keriK^LKi) = ker( yfLK 2 ) = {0} and the matrix KlLK 2 is 
invertible. 

Remark 2. a) Intuitively spoken, the degeneracy of the original Lagrangian is due to the lack of magnetic energy 
terms for the capacitors. With each independent constraint on the capacitor currents, one degree of freedom of 
the system can be removed. Hence, as many capacitors constraints are required to remove the capacitor current 
(cf. [10]). 
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b) In addition, for a RLC (resp. RCLV) non-degenerate circuit, the number of resistors (resp. and voltage sources) 
has to equal the number of independent constraints involving the currents through the resistor (resp. and voltage 
source) branches. 

Theorem 1 (Equivalence). The system (11) and the reduced system (17) are equivalent in the following sense: 

(i) Let (q,p,v) be a solution of the reduced system (17) and let q — K 2 q,v — K 2 v and (q,p) — ¥£(q,v). Then 
(q, v, p) is a solution to system (11) and it holds p — K^p. 

(ii) Let (q, v, p) be a solution to system (1 1) and q — K^q, v = K^v and let p = K^p with the well-defined pseudo- 
inverse ofK-i {with Ki,K 2 = I). Then (q,p,v) is a solution of the reduced system (17). 

Proof. (i) Assume (q,p, v) is a solution of (17). From the assumption p = ¥L(q, v) it follows p - ^£{q, v) = 0. 
With t M = K*£ it holds 

d-C M ^ d£,„_„ n (dqY d£,„_ „ T d£, 



'-(<!■ r) = ^(K 2 q,K 2 v) = (|?) ^(K 2 q,K 2 v) = Ki -=(>/. r). 
oq oq \dq) dq dq 



d£ M T dt 
Similarly, it holds — — (q, v) = K 7 —— (q, v) and thus it follows 

dv dv 

dt M T dl T 

P = -gz-(q, v) = K T 2 —{q, v) = K T 2 p. 

Together with (14), this gives 

T ■ d£ M M T ldt \ T [dL \ dt T 

With ker(Kl) = im(K), it follows 

^-p + f L€ im(K) = A Q (q) 
as can be seen from diagram (5). Furthermore, we have 

q = K 2 q = K 2 v = v. 

and since it holds v = K 2 v, it follows from diagram (5) K T v = 0. Both expressions are equivalent formulations 
oftheKCL. 



(ii) Now assume (q, p, v) is a solution of (11). With ker{K^) = im(K) and (14), it follows immediately 

p - K * p - K *\-dq +fL )--dq- +fL - 

d£, d£, A 

Furthermore, from q = v we get Ki q = Ki v which gives q-v. Finally, we have p = K\p — K\ 



dv dv 

Remark 3. We require the assumption (q, p) = ¥L(q, v) (the fulfillment of the Legendre transformation) in Theo- 
rem l(i) for the fulfillment of the relation (11c). A unique derivation of p directly from p is in general not possible 
from p = K^p as it is for q and v: Although there is a canonical projection K^ : T*Q — > T*M, there is no correspond- 
ing canonical embedding of T*M into T*Q (see also [30]). Assuming p - K^p instead of (q,p) = ¥L(q,v) in (i), we 
only get the relation K^(p - d-C/dv) = 0, and (1 lc) may not be fulfilled. 
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4. Discrete variational principle for electric circuits 



In this section, we derive a discrete variational principle that leads to a variational integrator for the circuit sys- 
tem. Since the solution of the reduced system (17) can be easily transformed to a solution of the full system (11) 
(Theorem 1), we restrict the discrete derivation to the reduced case. For the case of a degenerate reduced system, the 
choice of discretization is important to obtain a variational integrator that manages to bypass the difficulty of intrinsic 
degeneracy and thus is applicable for a simulation. In this section, three different discretizations are introduced that 
result in three different discrete variational schemes for which the solvability conditions are derived. 

For the discrete variational derivation, we introduce a discrete time grid At = {t k = kh \ k = 0, . . . , AO, Nh = T, 



where N is a positive integer and h the step size. We replace the charge q : [0, T] 
and the flux linkage p : [0, T] — > TIM by their discrete versions q d : {t k \ 



N 

k=0 



M, the current v : [0, T] -> T- q M 



M, v d 



{tkf 



k=0 



TqM and 



Pd '■ i f i)f = o — * T~M, where we view q k = qAkh), = Vd(kh) and p^ = pd(kh) as an approximation to q(kh), v(kh) and 
p(kh), respectively. 



4.1. Forward Euler 

We replace the reduced Lagrange-d'Alembert-Pontryagin principle with a discrete version 
Slhj] (£ M (Vk, h) + (Pk, q -^-^ - vty i+hj] f»(q k , v t , t k )6q k = 0, 



k=0 



(20) 



k=Q 



where in (20) the time derivative q(t) is approximated by the forward difference operator and the force evaluated at 
the left point. 

For discrete variations 6q k that vanish in the initial and final points as dqo = 5q^ = and discrete variations 6v k 
and 6p k this gives 



M 



dv 



i N-\ 

(qo,v Q )- p Q ,Sv Q \ + ^ 

I k=l 



d£ M 1 

-Q~-(Mk, h) - jiPk - Pk-i) + fiiqk, h, h), 6q k 



I c~ q k - g k -i - \^/^ M 



(qk,v k ) - p k ,6v k 



. x . qN-qN-\ . , „ 
+ \0PN-1, ; V N -i ) = 0. 



(21) 



This leads to the discrete reduced constrained Euler-Lagrange equations 



M 



dv 



'—(qo,vo) = po 



d£, M 1 

-^r(qk, Vk) - j(Pk ~ Pk-i) + fiiqk, h, h) = 
dq h 

qk - qk-i 



d£ 



M 



h 



= Vk-l 



dv 



<qk,Vk) = Pk 



k= l,...,N-l 



qN - qN-i 
h 



V N -\. 



(22a) 



(22b) 



(22c) 



For the Lagrangian defined in (13) and the Lagrangian forces defined in (14), this results in 



po = K 2 LK 2 v 

- Y*"' = K l i- CK 2qk - diag(R)K 2 v k + u s (t k )) ' 
qk - qk-i _ \k=l,...,N-l 



h 



V k -1 



K^LK 2 v k = p k 



qN - qN-i 



= Vjv-i. 



(23a) 



(23b) 



(23c) 
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This gives the following update rule: For given (qo, Vo), use (23a) to compute po- Then, use the iteration scheme 



r 
*-2 



/ 





k!lk 2 



to compute q\,.. . ,<?jv, v\,...,v N and pi,... ,pN- 



\ 


'qk 




i 


hi 


0^ 




'q k -i 




' N 




h 















Vjfc-l 


+ 





) 






.0 





I) 




,Pk-l, 




Ml, 



u s {t k ) for k = l,...,N 



(24) 



Proposition 2. System (24) z's uniquely solvable if the matrix K^(L + hdiag(R))K2 is regular. 



Proof. System (24) is uniquely solvable if the iteration matrix A 





KlLK 2 



0\ 
-I 



JiKlCK 2 hKldi&g(R)K 2 I ) 



has zero nullspace. 



For Az = for z = (q, v, p), it holds (i) q = 0, (ii) p = AT* ZJsT 2 v, (iii) /z^ CK 2 q + hK' 2 &mg(R)K 2 v + p = 0. Substituting 
(i) and (ii) in (iii) gives K[(L + hdiag(R))K 2 v = 0. Thus z = is the unique solution of Az = iff K%(L + hdmg(R))K 2 
has zero nullspace. 

4.2. Backward Euler 

Approximating the time derivative q(t) by the backward difference operator rather than by the forward difference 
operator as 

IN _ 

q k ~q k -i . \|| \ ,-Mir. :-. , ,;- _n (25) 



with discrete variations 8q k , that vanish in the initial and final points as Sqo = 8q^ = and discrete variations 5\>k and 
8p k yields 



+ 



(sp u ^ - n) + (^(*.*0 - ft.*) + g [(%. ^ - n) 

I d£ M 1 \ I d£ M 

l-Q~-(qk-\,Vk-i) ~ y^Pk- Pk-i) + fL(qk-uv k -ut k -\),Sq k -iJ + l-^-(qk,h) - Pk,Svk 



(26) 



= 0. 



This gives a slight, but in this case significant, modification for the Euler-Lagrange equations as 

q\ - qo 



h 



= v\ 



iquvi) = pi 



1 



—^riqk-uh-i) ~ j(Pk -Pk-i) + f¥(qk-uv k -utk-\) 
oq h 

qu - qk-\ 



Vk 



~~dv~ 



(qk,Vk) = Pk 



k = 2,...,N 



(27a) 
(27b) 

(27c) 



Note that in contrast to the variational scheme (22) consisting of an explicit update for the charges q and an implicit 
update for the fluxes p, we now get an implicit scheme for q and an explicit scheme for p. In particular, for the 
Lagrangian (13) and the forces (14), we obtain the following update rule: For given (qo,vo) compute po via po = 



KILK 2 vq. Then, use the iteration scheme 



7 


-hi 


N 


'qk 


1 







-/ 


h 




.0 





/, 


,Pk, 


V 








\ 


'qk-i 




' ^ 




h-i 


+ 







,pk-l) 




JiK T 2 , 



uM-i) fork=l,...,N (28) 



to compute q\,.. - ,qN, V\,...,v N and p\,...,p N . 
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Proposition 3. System (28) is uniquely solvable if the matrix K^LK^ is regular. 



(I -hi 



Proof. System (28) is uniquely solvable if the iteration matrix A 



0\ 



KlLK 2 -I 











I) 



has zero nullspace. For Az = 



with z = (q,v,p), it holds (i) q = hv, (ii) p = K^LK^v, (iii) p = 0. Thus z = is the unique solution of Az = iff 
K^LK 2 has zero nullspace. 

Proposition 3 says that whenever the KCL cancels the degeneracy of the system, the backward Euler scheme is ap- 
plicable, whereas the forward Euler scheme is applicable to a wider class of circuit systems (cf. Proposition 2) for h 
sufficiently large. The resulting variational Euler schemes (22) and (27) consisting of a combination of implicit and 
explicit updates are first order variational integrators. The construction of higher order implicit schemes (e.g. varia- 
tional partitioned Runge-Kutta (VPRK) methods along the lines of [5]) allows the simulation of arbitrary circuits. As 
an example, we present in the following a variational integrator based on the implicit midpoint rule. 



4.3. Implicit Midpoint Rule 

We introduce internal stages Q k ,P k ,V k , k = l,...,N - 1 that are given on a second time grid At = {t^ = 



(k + ±)h\k 



Pd : {T k 



0, . . . ,N - 1} and define the internal stage vectors : {r k ) 



k=o 



M, V, 



™*=0 



T- q M and 



,N-1 
<k=0 



TIM to be V k = v(tk + jh), Q k 



qk + \hv k ,P k 



ffi 



(Qk, V k ). The approximations at the nodes are 



then determined by the internal stages via q k+ \ = q k + hV k and p k+ \ - p k + h^-(Q k , V k ). 

Taking variations 6q k ,6Q k ,6p k ,6P k ,6V k for the following discrete Lagrange-d'Alembert-Pontryagin principle 
with 6qN = but free Sqo and initial value q Q 



N-i 



s\hJ](£ M (Q k ,V k ) + (p k , - M + (p k+l , - V k )\ + (po, qo ~ q°) 



k=0 



h 

N-l 



+hYjfL(Qk,V k ,T k )8Q k = Q 



k=0 



(29) 



gives 



+ (SP k , 



N-l 

2 

k=Q 

Qk 



-=H6*, %) + -f + f?{Qk, Vk, n), 5Q k \ + i-j^iQk, v k ) ■ 



^Pk~ Pk+l,SV k 



qk 



1 V k ) + (sp. 



h 2 
The Euler-Lagrange equations are 



k+l, 



qk+ i - qk 

h 



M- 



■ Pk+l + Pk 



h 



,Sq t 



+ (dp ,qo-q°) = 0. 



(30) 



-=r-(&, V k ) + ^+ f?(Q k , V k , T k ) 

aq h 

dt M - - 1 - 
"3=-(G*. V k) ~ 2 Pk - 

Qk-qk 1 f7 
—h—'2 Vk 

qk+i - qk 



h 



- V k = 

-Pk-Pk + i+Pk = 0, k = 0,...,N-l 
qo-q° = 0. 



(31a) 
(31b) 
(31c) 

(3 Id) 

(31e) 
(3 If) 
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Eliminating P k by equation (31e) together with V k = v k+ i, Qk = gt+ |* +1 (which follows from (31c) and (3 Id)) and 
Tk = '^^r 1 = t k+ i_ leads to the iteration scheme 

Pih-i =Pt + A-^-[ 2 + \ ^ .v t+ J (32a) 

ft+i = ft + hv k+ x_ (32b) 
pk+pk+i d£ M (q k +q k+ i _ \ 



Remark 4. The integrator (32) is equivalent to a Runge-Kutta scheme with coefficients a - \,b - \,c - \ (implicit 
midpoint rule integrator) applied to the corresponding Hamiltonian system. 

For the circuit case with Lagrangian (13) and forces (14), we start with given (qo,po) to solve iteratively for 
(q k +i,v k+ i, pk+i), k = 0, . . . , N - 1 for given u s (t) using the scheme 



-hi 

k!lk 2 



\\hKlCK 2 hK^diag(R)K 2 
fork = 0,...,N-l. 



I 



\ 


'qk+\ 










/ 


,Pk+l) 
















' N 






+ 







I Pk ) 







(33) 



Remark 5. The discrete current v k+ i, playing the role of the algebraic variable in the continuous setting, is only ap- 
proximated between two discrete time nodes t k and t k+ \. Also, note that v k _i is not explicitly used for the computation 
of (q k +i, v k+ i ,p k +i) (which corresponds to a zero column in the matrix of the right hand side of (33)). This means, that 
the computation of the magnitudes at time point t k+1 depends only on the discrete magnitudes within the time interval 
[tk, h+i], which is characteristic for a one-step scheme. In particular, v_ i (k = 0) is a pseudo-variable that is not used. 

Proposition 4. System (33) is uniquely solvable if the matrix Kl(2L + hdiag(R) + lh 2 C)K2 is regular. 



Proof. System (33) is uniquely solvable if the iteration matrix A = 



I 


{\hKlCK 2 



-hi 

K^LK 2 
hKldiag(R)K 2 



~\1 
I ) 



has zero 



nullspace. For Az = with z = (q, v, p), it holds (i) q = hv, (ii) p = 2K^LK 2 v, (iii) \hK\CK 2 q + hK\&\&g(R)K 2 v +p = 



0. Substituting (i) and (ii) in (iii) gives K^(2L + hdiag(R) - 
if K[{2L + hdmg(R) + \h 2 C)K 2 has zero nullspace. 



\h 2 C)K 2 v = Thus, z — is the unique solution of Az = 



Note that for linear circuits, the condition given in Proposition 4 is fulfilled for h sufficiently large, if the continuous 
system (18) has a unique solution. 

Remark 6 (Condition numbers). From Proposition 1, we see that the continuous reduced system (17) may be still 
degenerate due to the intrinsic degeneracy of the circuit topology and configuration. On the discrete side, both the 
forward Euler and the midpoint integrator show some regularization property: by perturbing K^LK 2 in magnitude 
proportional to the time step size h, both integrators can render the degenerate continuous reduced system (17) into 
regular discrete systems (24) and (33), respectively. However, this regularization comes at the price of possible large 
condition numbers as explained in the following. The iteration matrices A of the different schemes can be written 

7 



as A = A + hE with A = 



K{ LK 2 -I 
/ 



and E given by the respective iteration scheme. If the reduced 



system is regular (i.e. K\LK 2 is non singular), Aq is non singular with positive constant condition number k(Aq) and 
k(A) approaches a positive constant when the step size h goes to zero (e.g. one can compute the singular values of 
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the perturbed matrix A + hE using arguments from perturbation theory). In this case all iteration schemes are well 
conditioned independent of the step size h. However, if the reduced system is degenerate, i.e. K^LK^ is singular, 
also Aq is singular and the condition numbers of the forward Euler and the midpoint scheme grow reciprocally to the 
time step size h, i.e. k(A) « 0(1 1 h) for small h. When the circuit topology is fixed, the circuit's physical parameters 
are constants and the time step h is also fixed, preconditioner can be precomputed and applied to the systems (24) 
and (33) to improve their numerical stabilities. Since this work mainly focus on the theoretical aspects of variational 
integrators, we leave this stability issues for future work and assume thereafter no such issues in the subsequent 
discussion. 

5. Structure-preserving properties 

In this section, we summarize the main structure-preserving properties of variational integrators (see e.g. [30]) and 
their interpretation for the case of electric circuits. 

5.1. Symplecticity and preservation of momentum maps induced by symmetries 

Symplecticity. The flow on T*Q of the Euler-Lagrange equations preserves the canonical symplectic form of the 
Hamiltonian system. Variational integrators are symplectic, i.e. the same property holds for the discrete flow of the 
discrete Euler-Lagrange equations; the canonical symplectic form is exactly preserved for the discrete solution. Using 
techniques from backward error analysis (see e.g. [15]), it can be shown that symplectic integrators also have good 
energy properties, i.e. for long-time integrations, there is no artificial energy growth or decay due to numerical errors. 
This can also be observed for our circuit examples in Section 7. In the case of linear LC circuits that involve a 
quadratic potential, a second order variational integrator, e.g. the midpoint variational integrator as derived in Section 
4.3, even exactly preserves the energy (magnetic plus electric energy). 

Preservation of momentum maps induced by symmetries. Noether's theorem states that momentum maps that are 
induced by symmetries in the system are preserved. More precisely, let G be a Lie group acting on Q by (D : G X Q — > 
Q. We write <b g := 0(g, ■). The tangent lift of this action O rG : G x TQ -» TQ is given by <& T g Q (v q ) = T(O g ) ■ v q with 
v q e TQ. The action is associated with a corresponding momentum map J : TQ — > q*, where g* is the dual of the Lie 
algebra q of G. The momentum map is defined by 

{J(q, v), = | ^ , Z Q {q) J = (p, i; Q (q)) V£ € Q, 

where £ e is the infinitesimal generator of the action on Q, i.e. := j t \ _ Q Q>(exp(t£;),q) and exp : g — > G is the 

exponential function. Having a holonomic system described by a Lagrangian £, and a holonomic constraint h(q) = 0, 
then this system has a symmetry if the Lagrangian and the constraint are both invariant under the (lift of the) group 
action, i.e. X ° = X and h o (£> g (q) = for all g e G. Noether's theorem states, that if the system has a symmetry, 
the corresponding momentum map is preserved. In presence of external forces, this statement is still true, if the force 
is orthogonal to the group action. The discrete version of Noether's theorem [30] states, that if the discrete Lagrangian 
has a symmetry, the corresponding momentum map is still preserved. The variational integrator based on this discrete 
Lagrangian is thus exactly momentum-preserving. For constrained and forced systems, the preservation still holds 
with the additional invariance and orthogonality conditions on constraints and forces analogous to the continuous 
case. 

Considering an electric circuit, we are faced with a constrained distribution given by the KCL and external forces 
due to resistors and voltage sources. The KCL are formulated on the tangent space; however, since these are linear, 
they are integrable resulting in KCL on the configuration space. Thus, in the following, we are able to apply the theory 
of holonomic systems to derive a preserved quantity for an eletrical circuit under some topology assumptions of the 
underlying graph. 

Proposition 5 (Invariance of Lagrangian). The Lagrangian (6) of the unreduced system is invariant under the trans- 
lation of q L . 
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' v L - 




' qi + g " 




' qi + gi ' 


vc 


L 


vc 


1 


qc 
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qc 


vr 




Vr 


~ 2 


qR 




qR 


, Vy ) 




, Vy ) 




v qv ) 




v qv ) 



Proof. Consider the group G = R" L with group element g e G. Let <J> : G x Q -> Q be the action of G 
defined as <£(g,#) = (q L + g,qc) for each g e G with tangent lift O rG : G x TQ -> TQ, <D rG (g, v)) = 
(<?l + g, <?c, <?v, v L , v c , v R , v v ). Then, it holds 



£o® T g Q(q,v)= - 

= l - V T Lv- l -q T Cq = Z{q,v), 
since C = diag . . . , with the first n L diagonal elements being zero. 

Assumption 1 (Topology assumption). For every node j, j = l,...,m in the circuit (except ground), the same 
amount of inductor branches connect inward and outward to node j. 

In particular, Assumption 1 implies, that the sum of each row of K T L is zero, i.e. = for i = 1, . . . ,m. 

Proposition 6 (Invariance of distribution). Under Assumption 1, the KCL on configuration level are invariant under 
equal translation of q L . 

Proof. The group element g € G describing an equal translation of all components of qL can be expressed as g = al 
with a e R and 1 being a vector in R" L with each component 1. It follows 



( qi + g 
qc 
qR 

qv 



= K T q + K T L g = K T q + K[la = K T q, 



K T o $> g (q) = K T 
since the sum of each row of K T L is zero. 

Proposition 7 (Orthogonality of external force). The external force f L (8) is orthogonal to the action of the group 
G = W L being translations of q L . 

Proof. Let £ e g = R" L be an element of the Lie algebra. For the group action <b g (q) = (qL + g, qc, qR, qv), the 
infinitesimal generator can be calculated as 



?fl(9) 



dt 



<*Wf(<7) = "77 

(=o at 



r=0 



(q L + exp f£ 4c, q R , q v ) = (£, 0, 0, 0) 



It thus holds, 



</l,?q(9)> = <-diag(fl)v + diag(£) M , £>(<?)> = 0, 
since diag(/?) and diag(£) have zero entries in the first lines and columns. 

Theorem 2 (Preservation of flux). Under Assumption 1, the sum of all inductor fluxes in the electric circuit de- 
scribed by the Lagrangian (6), the external forces (8), and the KCL is preserved. 

Proof. From Proposition 5, 6, and 7 we know that the Lagrangian and the KCL are invariant under the group action 
Qgiq) = (qL + g, qc, qR, qv) with g e G - R nL and the external force f L is orthogonal to this group action. It follows 
with Noether's theorem, that the induced momentum map is preserved by the flow of the system. For the momentum 
map, we calculate 



{J(q,v),0 



dv 



,ta(q) 



Thus, the preserved momentum map is J(q, v) = p„ Lj , i.e. the sum of the fluxes of all inductors in the circuit. 
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Remark 7 (Proof based on Euler-Lagrange equations). An alternative proof can be derived based on the Euler- 
Lagrange equations in the following way. From (12a), it holds 



p L = K L A. 



For the time derivative of the sum of all inductors, it follows 



i Z p' = Z a = Z Z< = z ^ z^ = °- 

i=l i=l i=l 7=1 j=l i=l 




since with Assumption 1, it holds = Z^C^)^ = (^l)j'; f° r i = 1, . . . , m. Thus, P; is preserved. 

Remark 8 (Momentum map for reduced system). Also, for the reduced system described by the Lagrangian (13), 
the same momentum map can be computed by considering the group action ®g(#) = q + g with the group element 



Lemma 3 (Preserved momentum map). For any linear circuit described by the Lagrangian (6), the external forces 
(8), and the KCL, the momentum map defined by rf 'j£ with r\ e ker(K[) is preserved. 

Proof. Using the Euler-Lagrange equations, we see immediately 



since rf e ker(Kj) ± im(K L ) 3 K L A and thus, if |£ = const . 

The discrete Lagrangian system, including constraints and forces introduced in Section 4, inherits the same sym- 
metry and orthogonality property as the continuous system. Due to the discrete Noether theorem, the resulting 
variational integrators exactly preserves the sum of inductor fluxes under Assumption 1 (compare Section 7.4 for 
a numerical example). 

5.2. Frequency spectrum 

As can be observed in numerical examples (see Section 7), the frequency spectrum of the discrete solutions is 
much better preserved using variational integrators than other integrators. 

We want to analytically demonstrate this phenomenon by means of a simple harmonic Id oscillator. Assume the 
curves (q(t), p(t)) on [0, T] describe the oscillatory behavior of the system. Consider the discrete solution {(qk, Pk)} k=Q 
defined on the discrete time grid {tk) n k=i) with to = 0, = T and h = t k +\ - t k that is obtained from the one-step update 
scheme (qk+\,Pk+i) T = A{q k ,p k ) T , k = 0,...,N— 1, where A e R 2 ' 2 depends on the constant time step h. We assume 
that the discrete solution {(qk, Pk)} k=0 converges to the solution (q(t),p(t)) for decreasing h. Since this solution is 
oscillating and due to the convergence of the scheme, at least one eigenvalue A\ of A has to be complex (with nonzero 
imaginary part) for a small enough time step h. Since A e R 2 ' 2 the second eigenvalue A2 has to be complex conjugate 
to the first one. Thus, the corresponding eigenvectors are linearly independent and A is diagonalizable as A — QVQ 1 
with V = diag(/li,/L;). With the coordinate transformation (xk,yk) T = Q~ 1 (qk,Pk) T it holds (xk+uyk+i) 7 = V(xk,yk) T , 
i.e. x k+ \ = A\x k andyt +1 = A 2 y k - 

We demonstrate the preservation of the frequency spectrum for the Id oscillator in two steps: (i) We show that 
for a convergent scheme the update matrix A has two eigenvalues both of norm 1 if and only if the update scheme is 
symplectic. (ii) We show that methods defined by matrices with norm 1 eigenvalues preserve the frequency spectrum 
defined on different time spans. 

(i) "<=": Assume the scheme defined by A is symplectic, then det(A) = 1 (see e.g. [29]). It follows with A\ complex 
conjugate to A 2 (A 2 = A\): 1 = det(0 • det(V) • det^ 1 ) = Ai ■ A 2 = \Ai\ 2 = \A 2 \ 2 and thus = 1, i = 1,2. 
"=>": Assume A has two complex conjugate eigenvalues A\ = A* 2 with \A\\ = \A 2 \ = 1, i.e. we write A\ = e' 8 
and A 2 = e~' e with 6 eR and V = diag(e ,e , e~' 6 ). Note that 8 depends on the constant time step h that is used 



g g G c R"- m defined asg = K : 



- + ( 8 
2 \ 



with being the well-defined pseudo-inverse of K 2 . 



d T dX. 
dt^ dv 



= rfp = T] T K L A = 0, 
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for the discretization. Let J = I ^ 1. I be the canonical symplectic form and introduce the non-canonical 



-1 

symplectic form J = Q T JQ. We show that V preserves J, and therefore A preserves J, i.e. A is symplectic. 
Since J is skew-symmetric with zero diagonal, J is of the form ^ ^ q j w i tn A e R. It follows 



<r ifl j y —A j\ <r ie / \ -e-'V e A / \ -A 

(ii) Suppose that the discrete values x\, xi, . . . , x^ determined by the update scheme A are known, and admit the 
following discrete inverse Fourier transformation 

1 £ . (2m 



n-\ * 



foi , fe = l,...,N. 

N 



Consider a sequence of discrete points {Xk}^ =l that is shifted by one time step such that = x^+i = X\X k , k = 
l,...,N, i.e. \Xk}^ =1 approximates the solution on a later time interval than {xk}^ =v This admits the following 
discrete inverse Fourier transformation 



1 N l2ni \ 

X k = — ^/iix„exp(— knV k= l,...,N, 

n= 1 * ' 



i.e., X„ = X\X n . By the definition of the frequency spectrum, it holds X*X n = x* n X\A\X n = x* n \A\\ 2 x n = x* n x n , 
where the last equality relies on the symplecticity. Shifting the discrete solution arbitrary times, we see, that 
the spectrum will be preserved using different time intervals for the frequency analysis. This means that, in 
particular for long-time integration, a frequency analysis on a later time interval yields the same results as 
on an earlier time interval, which we denote by preservation of the frequency spectrum. The analysis for y 
follows analogously, and with the linear transformation Q the same holds for q and p. On the other hand, if 
\Ajj\ + 1, i, j = 1,2 (such as for non-symplectic or non-convergent methods), the frequency spectrum will either 
shrink or grow unbounded. 

Although the analysis was only performed for the simple case of a Id harmonic oscillator (in particular statement 
(i) is restricted to this case), we believe that for higher-dimensional systems, a similar statement as in (ii) can also be 
shown, which is left for future work. 

Relation to numerical results. For the numerical computations in Section 7, we perform a frequency analysis in 
two different ways: Firstly, we calculate the frequency spectrum on different time subintervals of the overall time 
integration interval [0, T]. This is directly connected to the analytical result of frequency preservation, i.e. we can 
observe that the spectrum is independent on the specific time interval using a symplectic integrator; however, using 
a non-symplectic method, the spectrum is damped calculated on a later time interval. Secondly, we use a fixed time 
interval [0, T] for the frequency analysis, but use different time steps resulting in different iteration matrices A. As we 
saw for the harmonic oscillator using a symplectic method, the magnitude of the two eigenvalues is independent on 
the time step h where it might de- or increase for increasing h for a non-symplectic method (e.g. using the explicit 
Euler method, the absolute value of the eigenvalues is 1 + 0(h 2 ) and the frequency spectrum would grow for larger h). 



6. Noisy circuits 

In this section, we extend the constructed variational integrator to the simulation of noisy electric circuits, for 
which noise is added to each edge of the circuit. Following the description in [6], in the stochastic setting, the 
constrained stochastic variational principle is 

6 f £(q(t\ v(f)) + (p{t), q(t) - v(0> dt + f f L (q(t), v(f), t) ■ Sq{i) dt + f 6q(t) ■ (2 o dW t ) = (34) 
Jo Jo Jo 
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with constrained variations 6q e A Q (q), where 2 is a n x n matrix, usually constant and diagonal, indicating the 
amplitude of noise at each edge, W, is a n-dimensional Brownian motion, and the last stochastic integral is in the 
sense of Stratonovich. This principle leads to the constrained stochastic differential equation 

|-^/ i+ Eo^ £ A»( ? ) (35a) 

dq = vdt (35b) 

8 f-P = (35c) 
ov 

K T v = 0, (35d) 

where by (35a) we mean that it holds f Y^-dt - dp + f L dt + £ o dWA = f X(q)dt for a vector field X{q) e A ( i (q) 

Jo \dq ) Jo J 

for any T. Correspondingly, the reduced stochastic variational principle reads 

£ M (q(t),v(t)) + (p(t),Ht)-m) dt+ fM( q (t),v(t),t)-Sq(t)dt+ 8q(t) ■ (^[Z o dW,) = 0. (36) 
o Jo Jo 

This results in the reduced stochastic Euler-Lagrange equations 

q_t m 

-j^rdt -dp + f^dt + K T {L o dW, = (37a) 

dq = vdt (37b) 

d£ M 



dv 



p = 0. (37c) 



To derive the discrete equations with noise, the Stratonovich integral is approximated by a discrete version. For 
simplicity, we present the equations for the forward Euler iteration scheme only. On the interval [tk, f/t+i] the integral 
j t 4+1 Sq(t) ■ (K^Y. o dW t ) is approximated by the discrete expression 8q k ■ (K^r)B k with B k ~ JV(0, h), k = 0, . . . , N - 1 
(see also [6]). In this way, we obtain the following reduced stochastic discrete variational principle 

(N-\ _ ^1 N-l N-l 

h ^ U M (q k , v*) + I p k , qk+l ~ qk - h)) h, t k )6q k +^hj] K ^ k ' ^ = °' (38) 

*=o v x h '' J k=0 k=0 

where for each k = 0, . . . , Af - 1, & is a n-dimensional vector with entries being independent standard normal random 
variables. The discrete reduced stochastic Euler-Lagrange equations that give the symplectic forward Euler iteration 
scheme is then given by 

d£. M 1 

-Q~~^k, h) - jiPk - Pk-i) + fiiVk, v k , t k ) + -^K^ k = 

q k ~ qk-i 

h = Vjt_1 

d£ M 

-^—(qk,vk) = p k , k=l,...,N. 
ov 

Different symplectic variational schemes (e.g. backward Euler or midpoint scheme) can be derived in the same way as 
in Section 4 with an appropriate discretization for the Stratonovich integral. In [6], it is shown that the stochastic flow 
of a stochastic mechanical system on T*Q preserves the canonical symplectic form almost surely (i.e., with probability 
one with respect to the noise). Furthermore, an extension of Noether's theorem says that in presence of symmetries of 
the Lagrangian, the corresponding momentum map is preserved almost surely. 

7. Examples 

In the following section, we demonstrate the variational integration scheme by means of simple circuit examples. 
The numerical results are compared with solutions resulting from standard modeling and simulation techniques from 
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circuit theory. In particular, we compare the variational integrator results based on Lagrangian models with solutions 
obtained with a Runge-Kutta scheme of fourth order as well with solutions obtained by applying Backward Differen- 
tiation Formula (BDF) methods to models derived using the Modified Nodal Analysis (MNA). For all methods, we 
use a constant step size h. 

For all examples, we use the convention, that the charge vector q e W is ordered as q — (qL,qc,lR,<lv), an d 
correspondingly the current, voltage, and linkage flux vectors as well as the Kirchhoff Constraint and the Fundamental 
Loop matrix. 



7.1. Short introduction to MNA 

In most circuit simulators, the Modified Nodal Analysis (MNA) is used to assemble the system of equations. In 
the following, we present the standard modified nodal analysis. We follow the description in [40]. A more detailed 
description can be found, for example, in [40, 14, 1]. The MNA consists of three steps: 1. Apply the Kirchhoff 
Constraint Law to every node except the ground. 2. Insert the representation for the branch current of resistors, 
capacitors and current sources. 3. Add the representation for inductors and voltage sources explicitly to the system. 

The combination of Kirchhoff's laws and the characteristic equations of the different elements yields the system 
of differential and algebraic equations 

K^qc(Kcuj) + Klg(K R u,t) + Klv L + Klv v +Kjv,(K U ,qc(Kcu,t),v L ,v v ,t) = (39a) 

p L {v L ,t)-K L u = (39b) 
Uy{Ku, qc(K c u, t), v L , v v , t) - Kyu = (39c) 

with node voltages u, branch currents through voltage and flux controlled elements v v and v L , voltage dependent 
charges through capacitors qc and current dependent fluxes through inductors pi, voltage dependent conductance g, 
and controlled current and voltage sources V/ and u v . System (39) can be rewritten in compact form as 



A[d(x(t),t)]' +b(x(t),t) = 



(40) 



with 





u 




\Kl 





X = 


vl 


, A = 





/ 




. v v _ 











d(x, t) 



q c (K c u, t) 
Pl(vl, t) 



and the obvious definition of b. The prime [d(x, t)\ = 4 [d(x(t), t)] denotes differentiation with respect to time. Since 
the matrix A ad< £^ might be singular, equation (40) is not an ordinary differential equation but of differential algebraic 
type. For a detailed description of the properties of these equations, we refer to e.g. [20]. 

The standard approach to numerically solve the system of equations (40) is to apply implicit multistep methods 
for the time discretization, in particular lower-order BDF schemes or the trapezoidal rule. For a detailed description 
of these methods, we refer to e.g. [40, 16]. The advantage of BDF methods is the low computational cost compared 
to implicit Runge-Kutta methods. However, they may have bad stability properties (cf. [40]) and, in particular, they 
are not symplectic. 



7.2. RLC circuit 

Consider the graph consisting of four boundary edges and two diagonal edges of a square (see Figure 2). On each 
edge of this graph, we have a pair of capacitor (with capacitance C, : = 1, i = 1, . . . , 6) and inductor (with inductance 
Li = I, i = 1,...,5) except on one edge. 2 On this edge, there is only one capacitor which leaves a degenerate 
Lagrangian. The corresponding planar graph consists of n = 6 branches and m + 1 =4 nodes, thus we have I = 3 
meshes. 

The matrix Kirchhoff Constraint matrix K e M" m and the Fundamental Loop matrix K2 e M"'"~ m are (with the 
fourth node assumed to be grounded) 



2 The values for inductance and capacitance are just chosen for demonstration purpose of the variational integrator. For real circuits, the values 
are typically of different order, resulting in a dynamical behavior on a different scale as in our numerical examples. 



20 



4 



Figure 2: Graph representation of a RLC circuit. 
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(41) 



The matrix K^LK^ is non-singular with L = diag(Li, . . . ,L$,0); thus, the degeneracy of the system is eliminated 
by the constraints on the system, and all three variational integrators derived in Section 4 can be applied. 

In Figure 3 a) the oscillating behavior of the current on the first branch is shown (the currents on the other branches 
behave in a similar way). For the energy behavior of the LC circuit, we compare the exact solution with solutions 
obtained with the three different variational integrators, a Runge-Kutta method of fourth order, and a BDF method 
of second order. Since no resistor or voltage source is involved, this energy should be preserved. For the variational 
integrator based on the midpoint rule (VI), the energy is exactly preserved since the electric potential is only quadratic. 
For this relatively short integration time span, we observe that also the solution with the Runge-Kutta scheme (RK4) 
preserves the energy (the red, blue and black lines in Figure 3 b) lie on top of each other). Using the forward (VI 
EFD, magenta) or backward Euler (VI EBD, green) variational integrator the energy oscillates around its real value, 
however, no dissipation or artificial growth of the energy occurs in contrast to the solution obtained by a BDF method 
(MNA BDF): Here, the energy rapidly decreases using a second order BDF method (see the cyan line in Figure 3 
b)). These results are based on a step size of h = 0.1. Increasing the step size to h — 0.4, a phase shifting of 
the currents computed with variational integrators is observed (which is a typical behavior observed for variational 
integrators) in contrast to solutions obtained by the Runge-Kutta scheme (see Figure 4 a)). 3 However, considering the 
energy behavior shown in Figure 4 b), energy is dissipating for the Runge-Kutta solution, whereas for the variational 
integrators it (its median, respectively) is preserved. The performance of the BDF method is even worse: The energy 



3 Note that the Runge-Kutta method is of higher order (fourth order) than the variational integrators, and thus solution curves are of higher 



accuracy. 
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Figure 3: LC circuit (no resistors) with step size ft = 0.1. a) The oscillating behavior of the current on the first branch is shown, b) Comparison 
of the exact energy behavior (exact) and the numerical solution using the three different variational integrators, midpoint rule (VI), backward Euler 
(VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order (RK), and a BDF method of second order based on MNA (MNA 
BDF). The energy is (qualitatively) preserved for VI, VI EBD, VI EFD, and RK. The use of BDF leads to an artificial energy decay. 



is dissipating very fast. This can also be observed considering the current in Figure 4 a). The amplitude of the current 
oscillations is damped to almost zero after a certain integration time. 

We investigate this phenomenon even more by comparing the amplitude of the oscillating branch current and the 
corresponding spectrum in frequency domain for the different integration methods. The solution of the first branch 
current oscillates with two frequencies, oji = 1 and a>2 = V2. In Figure 5, the branch current and the frequency 
spectrum for the discrete solution computed with two different step sizes (h = 0.1,0.4) and different integrators are 
compared. For a bigger step size h, the amplitude of the oscillations and the spectrum of the higher frequency are 
artificially damped for the Runge-Kutta and the BDF method. However, using a variational integrator, the frequency 
is slightly shifted, but the spectrum is much better preserved. For h = 0.4, we compute the frequency spectrum for 
three different time intervals, |o, [|, ^-J, and rj. Corresponding to our analytical result regarding frequency 
preservation (see Section 5.2), we see in Figure 6 that the frequency spectrum is the same independent on the integra- 
tion time using a variational integrator. However, using the Runge-Kutta or BDF method it is damped choosing a time 
interval after a longer integration time. 

In each branch of the circuit, we now add a small resistor with resistance R, = 0.001, i = 1, ... ,6. Again, we 
compare the oscillating behavior and the energy behavior of the numerical solution obtained by different integrators. 
Due to the resistors, the energy in the system decays. The rate of energy decay is shown for the exact solution in 
Figure 7 b). The variational integrator solution respects this energy decay much better than the Runge-Kutta and BDF 
scheme. 



7.3. Oscillating LC circuit 

As second example, we consider the LC circuit given in Figure 8. It consists of two inductors with inductance 
L\ — 1 and Li — 1 and two capacitors with capacitance C\ — 1 and C2 = 10 and thus has n — 4 branches and m + 1=3 



nodes. The Kirchhoff Constraint matrix K e 
node assumed to be grounded) 



and the Fundamental Loop matrix K2 e 
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Figure 4: LC circuit (no resistors) with step size h = 0.4. Comparison of the exact solution (exact) and the numerical solution using the three 
different variational integrators, midpoint rule (VI), backward Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order 
(RK4), and a BDF method of second order based on MNA (MNA BDF). a) The use of variational integrators (VI, VI EBD, VI EFD) leads to a 
phase shifting in the numerical solution of the current. With the BDF method, the oscillations are damped out. b) The energy is (qualitatively) 
preserved for VI, VI EBD, and VI EFD. The use of RK4 and BDF leads to an artificial energy decay. 



With «c - 2 and Kc 



having full rank, we can 




-1 
-1 1 

follow from Proposition 1, that the reduced Lagrangian system is 
non-degenerate, i.e. all three variational integrators derived in Sec- 
tion 4 can be applied. 

In Figure 9, the oscillating behavior of the branch currents on 
the inductors (a)-b)) and of the branch charges on the capacitors 
(c)-d)) is depicted. For a long-time simulation, we compare the q 
exact energy behavior of the LC circuit with the energy behavior 
of the solution obtained with the three different variational integra- 
tors (VI, VI EBD, VI EFD), a Runge-Kutta method of fourth order 
(RK4), and a BDF method of second order based on MNA (BDF 
MNA). As for the previous example, independent on the step size 
h, the energy is exactly preserved using VI based on midpoint rule, 
whereas the solutions using the Euler VI oscillate around the real 
energy value without dissipation or artificial growth of the energy 
(see Figure 10). Using the BDF method, the energy rapidly de- 
creases for increasing time step h. The solution of the Runge-Kutta 
scheme shows a similar decreasing energy behavior for the step 

size h = 0.1,0.2,0.4. However, for a step size of h — 0.6, the energy seems to converge to constant value that is not 
the exact energy value, but slightly lower. This might be due to the fact that the amplitude of the first lower frequency 
is almost preserved and only the amplitude corresponding to the higher frequency is damped out for increasing inte- 
gration time, as shown in Figure 1 1 for different time spans and a step size of h — 0.6. The same property is reflected 
by the plots in the frequency domain in Figure 12 for h = 0.4, where for increasing integration time the spectrum of 
the high frequency (a>2 ~ 1.43) is damped for the Runge-Kutta and BDF method and preserved using a symplectic 
method. This phenomenon is also confirmed by the frequency spectrum plot in Figure 13 for different step sizes. As 
for the previous example, the spectrum corresponding to the higher frequency is damped out for higher step sizes h 
for the BDF and the Runge-Kutta method. 



Figure 8: Oscillating LC circuit. 



23 




500 1000 1500 2000 2500 3000 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 

time cu 



h = 0.l h = 0.l 




h = 0.4 h = 0.4 



Figure 5: Current of first branch of LC circuit (left) and corresponding frequency spectrum (right). Comparison of the exact solution (exact) and 
the numerical solution using the three different variational integrators, midpoint rule (VI), backward Euler (VI EBD), and forward Euler (VI EFD), 
a Runge-Kutta method of fourth order (RK4), and a BDF method of second order based on MNA (MNA BDF). For increasing step size h, the 
damping of the amplitude of the oscillations and of the higher frequency spectrum increases using RK4 and BDF. For VI, VI EBD, and VI EFD, 
the frequency is slightly shifted, but the spectrum and the current amplitude is much better preserved. 



7.4. LC transmission line 



To demonstrate the momentum map preservation 
properties of the variational integrator, we consider 
the LC transmission line consisting of a chain of in- 
ductors and capacitors as illustrated in Figure 14. 
The matrix Kirchhoff Constraint matrix K e W ,m 
and the Fundamental Loop matrix K2 e H"*" - " 1 are 
(with the third node assumed to be grounded), 
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Figure 14: LC transmission line. 
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Figure 6: Frequency spectrum of first branch current of LC circuit (no resistors) (h = 0.4) computed on time interval [0, 773], [773,2773], 
[27/3, T]. Comparison of the exact solution (exact) and the numerical solution using the three different variational integrators, midpoint rule (VI), 
backward Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order (RK4), and a BDF method of second order based 
on MNA (MNA BDF). Using RK4 and BDF the spectrum is damped for higher integration times and preserved using a variational integrator. 



With n c — 2 and K c — | ^ ^ j having full rank, we can follow from Proposition 1, that the reduced Lagrangian 

system is non-degenerate. For this circuit, the topology assumption 1 holds, i.e. all nodes except ground have exactly 
one branch connected and inward and outward to the node. Thus, from Theorem 2 it holds that the sum of the inductor 
fluxes pi t + pi, + is a conserved quantity that can be also seen from the variational simulation results depicted in 
Figure 15. 

7.5. Validation on the stochastic variational integrator 
Consider the stochastic differential equation 

dx = Axdt + tdW t , (43) 

where £ is a n-by-ra matrix, not necessarily full rank, x = (jci.Xz, • • • , x„) e M", A e W'" and W, is an m-dimensional 
Brownian motion (with independent components). A common way to verify the quality of a numerical solution is to 
consider statistical moments of the solution; in particular, we focus on the expectation and the variance, i.e., K(x(t)) 
and D(x(0) = E(x(f)) 2 - (E(x(t))) 2 . 

On the analytical side, by Ito's formula (see e.g. [35]) we have with B(t) = exp(Af) 

E(x(0) = B(t)x(0) (44a) 
D(x(t)) = f B(t)11 t B(t) t dr. (44b) 



The expectation and the variance can always be computed if A and £ are given. For big systems, however, such a 
mundane computation is quite complex. On the numerical side, we run an ensemble of simulations (of total number 
M), all starting from the same initial condition but for each simulation an independent set of noise (i.e., different is 
used. The ensemble are indicated by x l (t), x 2 (t), . . . , x M {t) where for any j, x^(t) = (Xj(f), x J 2 (t), . . ., x J n (t)) is a vector. 
We compute the empirical moments by 
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Figure 7: LCR circuit (with resistors) with step size h = 0.4. Comparison of the exact solution (exact) and the numerical solution using the three 
different variational integrators, midpoint rule (VI), backward Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order 
(RK4), and a BDF method of second order based on MNA (MNA BDF). a) The use of variational integrators (VI, VI EBD, VI EFD) leads to a 
phase shifting in the numerical solution of the current. With the BDF method, the oscillations are artificially damped out. b) The energy decay is 
much better preserved for VI, VI EBD, and VI EFD as for RK4 and BDF. 

The numerical method is validated if for large enough M the empirical moments (45) are close to the analytical ones 
(44). 

In our setting, we can rewrite the reduced stochastic Euler-Lagrange equations (37) in the form of (43) with 

x = (q,p) e M 2( "~ m) , £ = | q j e M 2 (»-'»», and the obvious definition of A e r2(»-'»X2(»-») with I e E"-" 

and K 2 e W'-"-'". The analytical variance matrix B{{q{t),p(t)) e ]g}(n-m),2(n-m) for the re( i U ced system can now be 
calculated using equation (44b). The corresponding variance matrix for the full system can then be calculated as 



H(q(t),p(t)) = l K * ^ )m(t),p(t))l K i ®r V 



v,2n,2n 



As a demonstration, we calculate the empirical and analytical moments for the circuit introduced in Section 7.3. For 
the experiments throughout this section, we defined S as 4-by-4 diagonal matrix with diagonal entries = 0.01, j = 
1, . . .4. The step size is h = 0.1, the integration time for each simulation is T — 30, and we start with the initial 
conditions qo = (1,0), Vo = (0,0), and po = (0,0). The empirical averages are calculated over an ensemble of 
M = 100000 independent simulations. 

The analytical variance of p Ll and p Ll , i.e., the fifth and sixth diagonal elements of the variance matrix in the full 
system, are plotted as functions of time (see Figure 16, red dotted line). Notice, that p^ and p^ 2 in our case are just the 
currents through inductor branch 1 and 2, the inductances are L\ — L 2 — 1 . The result using the stochastic variational 
integrator is also shown in Figure 16 (blue solid line). Both function shapes and ranges agree very well. In particular, 
all the little bumps in the variance that are subtly different are approximated correctly. This classical test serves as an 
evidence that the stochastic integration works fine. 



7.6. Multiscale integration with FLAVORS 

When the circuit exhibits behavior in two time scales, our integrators can be FLAVORized [39] to capture the slow 
time scale without resolving the fast time scale to greatly reduce integration time. We first give a brief description of 
FLAVORS {FLow Averaging integratORS). For more details, we refer to [39]. 

Consider an ordinary differential system on W 1 



1 



if = G(u e ) + -F(u e ) 



(46) 
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c) d) 
Figure 9: LC circuit, a) Current on inductor 1. b) Current on inductor 2. c) Charge on capacitor 1. d) Charge on capacitor 2. 

with e <K 1 . In the context of Lagrangian systems, we consider a multiscale Lagrangian as 

£(q, V) = \ v T Lv - V(q) - - U(q) (47) 
2 e 

where V(q) is denoted as "slow" and - U(q) denoted as "fast" potential. In the case of a linear circuit, the slow potential 
corresponds to the charge potential of a capacitor with high capacitance, whereas the fast potential corresponds to a 
capacitor with very small capacitance. For instance, consider the oscillating LC circuit in Section 7.3 with C\ — 1 and 
Ci — e (e = 1 in Section 7.3). The corresponding potential in (47) can be written as the sum of slow and fast potential 
as 

V(q)=\q 2 Cl , U{q)= l -q 2 Cj _. 

The smaller e (i.e. Ci) is, the wider the two time scales will be separated. 

FLAVORS are based on the averaging of the instantaneous flow of the differential equation (46) with hidden 

slow and fast variables. The way to FLAVORize any of our integrators for circuits is to break each timestep into a 

composition of two substeps. The first one uses a timestep of length r, and the second one uses a timestep 8 — t.t has 

to be small enough to resolve the stiffness, but 6—r does not. In fact, in the first substep one integrates the entire system 

with the original value of the stiffness 1/e, whereas in the second substep one integrates the system with the stiffness 

turned off, i.e., 1/e is temporarily set to 0. Of course, for the first substep of the next step, 1/e has to be restored to 

1 

its original big value again. To be more precise, FLAVOR is implemented using an arbitrary legacy integrator for 
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h = 0.1 





h = 0.2 




h = 0.4 



h = 0.6 



Figure 10: Energy of oscillating LC circuit: Solutions obtained with the midpoint variational integrator exactly preserve the energy (VI) independent 
on the step size h. The energy computed with the forward (VI EFD) and backward (VI EBD) Euler variational integrators oscillates around the real 
energy value without dissipation or artificial growth. For solutions obtained with a second order BDF method (BDF MNA), the energy dissipates 
due to numerical errors. The solution of the Runge-Kutta scheme (RK4) dissipate energy, but seems to converge to lower constant energy value. 



(46) in which the parameter - can be controlled. By switching on and off the stiff parameter, FLAVOR approximates 
the flow of (46) over a coarse time step H (resolving the slow time scale) by the flow 



H - z 



where r is a fine time step resolving the "fast" time scale (t <k e) and M is a positive integer corresponding to the 
number of "samples" used to average the flow (S = H/M). Since FLAVORS are obtained by flow composition, they 
inherit the structure-preserving properties (for instance, symplecticity and symmetries under a group action) of the 
legacy integrator for Hamiltonian systems. 

Theorem 1.4 in [39] guarantees the accuracy of FLAVORS for 6 <K ho, t <k e and (~J « 5 « ^, where ho 
is the stability limit of step length for the legacy integrator. Furthermore, if the hidden fast and slow variables are 
affine functions of the original variables (such as in our case of linear circuit), the condition relaxes to 6 <K ho, t <k e 
and 5 <K 1 . A intuitive interpretation of that theorem is, the numerically integrated slow variable will convergence 
strongly (as a function of time) to the solution to an averaged effective equation, and the fast variable will converge 
weakly to its local ergodic measure. 



We use FLAVORS to simulate the LC circuit with e 



io- 



O.le = 10~ 4 , H = 0.1 and M = 100. The 



charges and currents as functions of time are plotted in Figure 17. Notice that the slow components in the solution are 
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Figure 1 1 : Charge on capacitor 1 of LC circuit (h = 0.6). Comparison of the exact solution (exact) and the numerical solution using the three 
different variational integrators, midpoint rale (VI), backward Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order 
(RK4), and a BDF method of second order based on MNA (MNA BDF). For the variational integrators the amplitudes for low and high frequencies 
are preserved. Using the Runge-Kutta method, the amplitude corresponding to the high frequency is damped out for increasing integration time 
whereas the amplitude of the lower one seems to be preserved. Using the BDF method, the amplitudes of both frequencies are damped. 

captured strongly, but the fast components may have altered wave shapes: for instance, Figure 18 shows a zoomed-in 
investigation of the current through the second branch, which is a superposition of a slow global oscillation and a 
fast local oscillation; the slow one is obviously well-captured in the usual sense, and the fast one is captured in the 
less-commonly-used sense of averaging. 

8. Conclusions 

In this contribution, we presented a unified framework for the modeling and simulation of electric circuits. Start- 
ing with a geometric setting, we formulate a unified variational formulation for the modeling of electric circuits. 
Analogous to the formulation of mechanical systems, we define a degenerate Lagrangian on the space of branches 
consisting of electric and magnetic energy, dissipative and external forces that describe the influence of resistors and 
voltage sources as well as (non-)holonomic constraints given by the KCL of the circuit. The Lagrange-d'Alembert- 
Pontryagin principle is used to derive in a variational way the implicit Euler-Lagrange equations being differential- 
algebraic equations to describe the system's dynamics. A reduced version on the space of meshes is presented that is 
shown to be equivalent to the original system and for which under some topology assumptions the degeneracy of the 
Lagrangian is canceled. 

Based on the reduced version, a discrete variational approach is presented that provides different variational in- 
tegrators for the simulation of circuits. In particular, the generated integrators are symplectic, preserve momentum 
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on [0.T/3J 



on [T/3,2T/3] 



on [2T/3.T] 




Figure 12: Frequency spectrum of charge on capacitor 1 of LC circuit (h = 0.4) computed on time interval [0, 773], [7/3,27/3], [27/3,7]. 
Comparison of the exact solution (exact) and the numerical solution using the three different variational integrators, midpoint rule (VI), backward 
Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method of fourth order (RK4), and a BDF method of second order based on MNA 
(MNA BDF). Using RK4 and BDF the spectrum is damped for higher integration times and preserved using a variational integrator. 



maps in presence of symmetries and have good long-time energy bahvior. Furthermore, we observe that the spec- 
trum of high frequencies is especially better preserved compared to simulations using Runge-Kutta or BDF methods. 
Having the variational framework for the model and the simulation, extensions of the approach using already-existing 
types of different variational integrators can be easily accomplished. As an example, we presented the extension for 
the simulation of noisy circuits using stochastic variational integrator approaches as well as multiscale methods (in 
particular FLAVORS) for an efficient treatment of circuits with multiple time scales. 

In the future, we will extend the approach to the simulation and analysis of more complicated nonlinear and 
magnetic circuits that might include nonlinear inductors, capacitors, resistors and transistors. Since a variational for- 
mulation in terms of energies, forces and constraints is still valid for the nonlinear case, the presented integrators will 
be derived and applied in straight forward way. Furthermore, the inclusion of controlled sources allows for the con- 
sideration of optimal control problems for circuits for which techniques also based on a variational formulation can 
be easily applied (see e.g. [34]). The variational simulation of combined mechanical and electric (electro-mechanical) 
systems is the natural next step towards the development of a unified variational modeling and simulation method for 
mechatronic systems. Furthermore, at nanoscales, thermal noise and electromagnetic interactions become an essential 
component of the dynamic of electric circuits. We plan to investigate the coupling of variational integrators for circuits 
with multi-symplectic variational integrators for EM fields and continuum mechanics (see e.g. [23, 28, 37]) to produce 
a robust structure-preserving numerical integrator for Microelectro mechanical and Nanoelectromechanical systems. 
Recently, Mike Giles has developed a Multilevel Monte Carlo method for differential equations with stochastic forc- 
ings [12] that shows huge computation accelerations, 100 times in some cases. The extension of the current method to 
multilevel stochastic variational integrators is straight forward and may further accelerate computation dramatically, 
especially for multiscale problems, while preserving certain properties of the circuit network. 
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Figure 13: Frequency spectrum of charge on capacitor 1 of LC circuit. Comparison of the exact solution (exact) and the numerical solution using 
the three different variational integrators, midpoint rule (VI), backward Euler (VI EBD), and forward Euler (VI EFD), a Runge-Kutta method 
of fourth order (RK4), and a BDF method of second order based on MNA (MNA BDF). For increasing step size h, the damping of the higher 
frequency spectrum increases using RK4 and BDF. For VI, VI EBD, and VI EFD, the frequency is slightly shifted, but the spectrum is much better 
preserved. 
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Figure 16: Benchmark of variances as functions of time according to (44b) (red) and variances as functions of time computed numerically by 
averaging over an ensemble according to (45b) (blue), a) Dpi b) Dp2- 
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Figure 17: Simulations of a multiscale system: a) Benchmark solution computed with a variational integrator (h = 10 4 ) b) FLAVOR with 
T = 1(T 4 ,<S = 1(T 3 and e = 1(T 3 . 
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Figure 18: Simulations of a multiscale system: a) Benchmark solution computed with a variational integrator (h = 1CT 4 ) b) FLAVOR with 
T = 1(T 4 ,<S= l(T 3 and£= 10~ 3 . 



35 



